#drafted by Jeremy VanDerWal ( jjvanderwal@gmail.com ... www.jjvanderwal.com )
#GNU General Public License .. feel free to use / distribute ... no warranties

################################################################################
################################################################################
#get the command line arguements
args=(commandArgs(TRUE)); for(i in 1:length(args)) { eval(parse(text=args[[i]])) }
print(ls())
# spp='15107'
taxa = c("Amphibians", "Birds","Mammals","Reptiles","Plants") ;	tax = taxa[1]
work.dir = paste("/home/ctbccr/Ramona/SDM3/",tax,"/models/",spp,'/',sep='') ; setwd(work.dir)

###############################################################################
#prep work and checks
library(SDMTools) #load the library
files = list.files(paste('output/ascii/'),pattern='asc.gz')

threshold = read.csv('output/maxentResults.csv'); threshold = threshold[which(threshold$Species==spp),]
spp.AUC = threshold$Training.AUC[1] #get the AUC
threshold = threshold$Equate.entropy.of.thresholded.and.original.distributions.logistic.threshold[1] #extract the species threshold value
dir.create(paste('summary/images/',spp,'/potential/',sep=''),recursive=TRUE); dir.create(paste('summary/images/',spp,'/realized/',sep='')) #create the output directories

#extract ES, GCM, year information
ESs = GCMs = YEARs = current = NULL
for (ii in 1:length(files)) {
varname = gsub('.asc.gz','',files[ii])		#.gz
tt = strsplit(varname,'_')[[1]] #string split the file name
if (length(tt)==1) { current = tt[1] } else { ESs = c(ESs,tt[1]); GCMs = c(GCMs,tt[2]); YEARs = c(YEARs,tt[3]) } #Assign the split string to the proper object
}
ESs = unique(ESs); GCMs = unique(GCMs); YEARs = unique(YEARs) #keep only unique values

################################################################################
clipasc = read.asc.gz('bioregion.clip.asc.gz') #read in prepared realised distribution clip

base.asc = read.asc.gz("/home/ctbccr/Ramona/Climate_Data/base.asc.gz")

pos = as.data.frame(which(is.finite(base.asc),arr.ind=TRUE))
pos$lat = getXYcoords(base.asc)$y[pos$col] #extract the longitudes
pos$lon = getXYcoords(base.asc)$x[pos$row] #extract the latitudes

pos$clip = extract.data(cbind(pos$lon,pos$lat),clipasc) #extract the clip information
pos = na.omit(pos)
cellarea = grid.area(base.asc) #get the area of individual cells
pos$area = cellarea[cbind(pos$row,pos$col)] #append the cell area
rm(clipasc)#; rm(occur); rm(regions) #clean up memory

################################################################################
#create a loop to cycle through RCPs to avoid overloading the memory
#es=ESs[1]; YEAR=YEARs[1]
#
for (es in ESs) {
	for (YEAR in YEARs) {
		# bring the data into memory
		tfiles = files[c(1,grep(es,files))] ; tfiles = tfiles[c(1,grep(YEAR,tfiles))]
		pot.mat = matrix(0,nr=nrow(pos),nc=length(tfiles)) #create matrix to store information
		colnames(pot.mat) = gsub('.asc.gz','',tfiles) #add column names  .gz
		
		for (tfile in tfiles) { cat(tfile,'\n'); 
			pot.mat[,gsub('.asc.gz','',tfile)] = read.asc.gz(paste('output/ascii/',tfile,sep=''))[cbind(pos$row,pos$col)] } #read in all the projection
		save(pot.mat,file=paste('output/',spp,'.',es,'.',YEAR,'.potential.dist.mat.Rdata',sep='')) #save the potential matrix
		pot.mat[which(pot.mat<threshold)] = 0 # change anything < threshold to 0
		real.mat = pot.mat; real.mat = real.mat[,] * pos[,'clip'] #create the realized distributions using the clip

		################################################################################
		# get some summary stats
		sum.data = function() { #this is a function to summarize the distribution data
		outdata1 = data.frame(ES=rep(NA,length(vois)),GCM=NA,year=NA) #define the basic information
		
		for (ii in 1:length(vois)) {
			tt = strsplit(vois[ii],'_')[[1]]
			if (length(tt)==1) { outdata1[ii,1:3] = c(NA,NA,1990) } else { outdata1[ii,1:3] = c(tt[1],tt[2],tt[3]) }
		}
		outdata1$sum.suitability = colSums(distdata*pos[,'area'],na.rm=TRUE) #get the sum of the environmental suitability
		outdata1$prop.abund = outdata1$sum.suitability / outdata1$sum.suitability[outdata1$year==1990] #calculate proportionate change in abundance
		#calculate the Class-based statistics and Istat, then append to the columns
		cur.asc = base.asc; cur.asc[cbind(pos$row,pos$col)] = distdata[,which(colnames(distdata)=="1990")] #define the current data surface for estimating Istat
		for (ii in 1:nrow(outdata1)) {
			tasc = base.asc; tasc[cbind(pos$row,pos$col)] = distdata[,ii] #put the data back into a matrix
			Ival = Istat(cur.asc,tasc) #calculate the Istatistic
			tasc[which(tasc>0)] = 1 #now convert binary data
			CS = ClassStat(tasc,latlon=TRUE) #get the class stats
			if (1%in%CS$class) { CS = CS[which(CS$class==1),] } else { CS = CS[1,]; CS[1,] = NA } #only keep info on distriubtion... if no distriubtion, set everything to 0
			if (ii == 1) { cois = NULL; for (jj in c('Istat',names(CS)[-1])) {outdata1[jj] = NA; cois = c(cois,which(names(outdata1)==jj)) } } #if the first summary, create columns to store data and define the column numbers for this data
			outdata1[ii,cois] = c(Ival,CS[,-1])
			}
			return(outdata1)#return the output
		}
		vois = colnames(pot.mat)
		distdata = pot.mat
		outdata = data.frame(dist.type='potential',sum.data())#extract the output summary data
		distdata = real.mat
		outdata = rbind(outdata,data.frame(dist.type='realized',sum.data())) #get the realized summary data
		
		write.csv(outdata,paste('summary/',spp,'.',es,'.',YEAR,'.classstats.csv',sep=''),row.names=FALSE) #write out the data
	}
}








